Labyrinthic Granular Landscapes 
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We have numerically studied a model of granular landscape eroded by wind. We show the ap- 
pearance of labyrinthic patterns when the wind orientation turns by 90°. The occurence of such 
structures are discussed. Morever, we introduce the density Uk of "defects" as the dynamic param- 
eter governing the landscape evolution. A power law behavior of Uk is found as a function of time. 
In the case of wind variations, the exponent (drastically) shifts from 2 to 1. The presence of two 
asymptotic values of Uk implies the irreversibility of the labyrinthic formation process. 
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I. INTRODUCTION 

The ripple formation due to the wind blowing across a 
sand bed has recently received much attention in the 
statistical physics community [^J^. Indeed, the physical 
mechanisms involved are complex phenomena of granular 
transport. Experi ental works as well as natural obser- 
vations have underlined the primary role played by 
saltation in the emergence of ripples and dunes. Along 
this line, various models for ripple formation have been 
proposed in the past. Theoretical models which consider 
hopping and rolling grains have generally led to travel- 
ing ripple structures Simulations have also 
considered various additional effects like the screening of 
crests, grain reptation and the existence of a grain ejec- 
tion threshold. Among others, the Nishimori-Ouchi (NO) 
model is able to reproduce a wide variety of different 
eolian structures: transverse, barkhanic, star-like, etc... 
The main advantage of such numerical model is that all 
parameters can be easily tuned and the physical mecha- 
nisms can be deeply investigated. 

In 1991, Goossens performed an original experiment 
which is the following. In a wind tunnel, a 12 x 12 crn^ 
rough granular landscape (dust size « 32 /i to) is eroded 
by an air flow (v ~ 132 cm s~^). This created small rip- 
ples perpendicular to the wind direction. After 45 min- 
utes, the wind orientation was changed by 90° and this 
poduced diagonal structures instead of a new set of rip- 
ples perpendicular to the previous ones. The Goossens's 
experiment represents a good test for the NO model. In 
the present paper, we report simulations of this partic- 
ular kind of landscape. This allows us to discuss the 
dynamics of this unusual phenomenon. 



II. MODEL 

In the Nishimori-Ouchi model Q, two kinds of granu- 
lar transport processes are considered: (i) the saltation 
and, (ii) the potential energy relaxation. The temporal 
evolution equation of the height of sand h{x) at point of 
coordinate x reads 



^d^h{x,t) 
dx"^ 



(1) 



At the right hand side of this equation, the first term rep- 
resents the saltation process. Due to wind shear stress, 
grains are moved from a position ^{x) to a position x. 
The mean amplitude of the path length is given by the 
constant A. On the other side, the constant _D is a relax- 
ation coefficient. This second term takes in account the 
transport phenomena along the slopes of the surface, e.g. 
reptation and avalanches. 

The NO model has been implemented as follows. A 
two dimensional square lattice with periodic boundary 
conditions is considered. To each site i,j of the lattice is 
associated a real number hi^j which represents the height 
of the granular landscape at that position. Assume that 
the wind blows along the i-axis. At each discrete time 
t, a site i,j is randomly choosen and a quantity qi j of 
matter is displaced by saltation from this site towards 
the site i + which is incremented by the height qi j. 
Both quantities £ij and qij are determined by 



£ij — a(tanh V/iij- + 1) 
= /^(l + £ ^ tanh V/iij) 



(2) 



where a and (3 are dynamical constants and the parame- 
ter e is the minimum quantity of sand which is displaced 
by saltation. The mathematical form (tanh Vft.) of those 
relationships assumes that the local slope mainly con- 
trols the granular transport. The flux of sand extracted 
from the faces exposed to the wind is indeed smaller than 
that screened by crests. After the saltation process (||), 
a relaxation of the landscape is assumed (creeping and 
avalanches) before the next time step t + 1 takes place. 
The relaxation reads 
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h^j{t + l) = h^j{t) 

+ D{^ Enn hnn{t) + ^ Ennn hnnn{t) - Kj{t)) 

+ 1) = hi+U,j.o{t) 

+ J2nn ^nn{t) + J^nnn hnnn{t) " h,+i^ ^^j{t)) 

where the summations run over nearest neighbors (nn) 
and next nearest neighbors (nnn) of both sites i,j and 
i-\-£ij,j. This equation is the discrete counterpart of the 
Laplacian relaxation of Eq. . The process is repeated a 
large number of times. Typically, we stop the simulation 
after t = 2.5 x lO"^ steps on 201 x 201 lattices. We inten- 
tionaly choose a lattice size that is not commensurable 
with the mean saltation lenght, a. 

III. RESULTS AND DISCUSSION 

We have performed extensive simulations by varying 
all the parameters: a, (3, D and e. Modifying a changes 
the mean ripple wavelength (the distance between two 
sucessive crests) , while D affects the aspect ratio (ampli- 
tude) of the ripples. The values taken by (3 and e permit 
or not the appearence of ripples: typicaly f3 S [0.2,0.6] 
and e = 0.3. Note that we have normalized the time t 
by the duration of the simulation, involving t G [0, 1] in 
arbitrary units (a.u.). 

Figure |l| shows a typical result of our simulations for 
a = 2.5, /3 = 5, D = 0.4, e = 0.3. The granular land- 
scape is shown for 4 different stages of evolution. One 
observes on the top row the early formation of ripples 
perpendicular to the wind direction. On the bottom row, 
the wind direction has changed by 90° clockwise and a 
labyrinthic structure appears. This observation empha- 
sizes the impact of the initial topography on the orien- 
tation of the ripples crests. One should note that such 
labyrinthic structures are strikingly similar to Goossens' 
ones Q. 




FIG. 1. Four different stages of a granular landscape evolu- 
tion within the NO model. When the wind direction changes, 
a labyrinthic pattern appears. The simulation parameters are: 
a =2.5, 13=5, D=OA, e = 0.3. The lattice size is 101x101. 
(a) t=0.2 a.u. and wind direction is up, (b) t=0.48 a.u. and 
wind direction is up, (c) t=0.52 a.u. and wind direction is 
right, (d) t=l a.u. and wind direction is right. 

In order to quantify the effect(s) of wind variations, 
we have measured the maximum ripple amplitude A^ax 
for both constant and variable winds. This quantity was 
also experimentally measured in j^]. As the wind orien- 
tation is changed by 90°, no significant change of Amax is 
observed, similarily to experiments This means that 
a brutal change in the wind direction does not modify 
the net deposit of sediments on the crests. The compe- 
tition between transport and deposit phenomena is not 
deteriorated by the perturbation of the wind orientation. 
Actually, A^ax represents adequately this competition, 
but is not a relevant dynamical parameter in order to 
understand the formation of labyrinthic structures. 

Looking for details in Figure |l], one can observe that: 
(i) diagonal structures appear at the vicinity of "defects" 
of the primary landscape. Those "defects" are kinks and 
antikinks. A kink is a bifurcation of a crest, while an an- 
tikink is a termination of a crest, i.e. a bifurcation of the 
valleys. Moreover, kinks and antikinks are not indepen- 
dent at all: they are formed by pairs. Kinks and antikinks 
can be considered as nucleation centers for new ripples 
when the wind direction changes. One understands the 
formation of labyrinths as follows. Old ripples are pushed 
in the new wind direction. If their crests are perpendicu- 
lar to the new direction ripples are compressed. However, 
near a "defect" , the angle between the crest and the wind 
is smaller than 90°. A rotation of the crest is thus ini- 
tialized there. This leads to diagonal structures, (ii) 
The formation of a labyrinthic-like structure involves the 
growth of the number of "defects" . 

Let us consider the relevant parameter: the density nfc 
of kinks. This quantity is defined as the number of kinks 
present on the surface, divided by the area of the lattice. 
In order to measure the number of kinks present in the 
landscape, we proceeded as follows (see the illustration in 
Figure |^) . The surface is recorded in grayscale images at 
different stages of evolution. The darkness indicates the 
height of sand, i.e. crests are in black while valleys are 
in white. Images are then analysed using common tools 
of image analysis. First, a threshold is applied in order 
to get binary (black/white) images with crests in black. 
Then, a function reduces all crests to a skeleton through 
an iterative erosion technique. The last step concerns the 
countdown itself. The program browses the skeleton line 
by line. When a black point is met (a crest), the number 
of its black neighbors is counted. If this number is greater 
or equals to 3, the point is necessary a kink. One should 
note that this method can be applied to images of real 
experiments. 
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FIG. 2. Illustration of the procedure of kink counting in 
the case of constant wind (top row) and variable wind (bot- 
tom row). Grayscale images are created with crests in black 
(left). The binary (black/white) images are extracted and 
show the crests (midle). An iterative erosion technique leads 
to a skeleton (right). 

Figure ^ presents the temporal evolution of n/j in both 
cases: constant and variable wind. One can see that 
without any wind modification, decreases as a func- 
tion of time (black squares). When a wind change occurs 
at time T, the density of kinks suddenly increases. In 
Figure ^ different wind variations arc illustrated for dif- 
ferent times T . After the jump of observed at time T, 
the kink density continues to decreases slowly. 
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FIG. 3. The density of kinks as a function of time t. The 
simulation parameters are : a =2.5, /3=5, D=0.4, e = 0.3, the 
lattice size is 201x201. Different times T for wind orientation 
changes are illustrated. The curves are fits for both cases: 
constant wind (bottom curve) and variable wind (top curve). 



the exponent d captures the dynamics of the landscape. 
Indeed, a large value of d involves a fast decrease of the 
kink density; such a situation implies that the surface 
can be easily modified by the wind. On the other hand, 
a small value of d means that the patterns are less af- 
fected by the variation of the wind orientation. Looking 
for details in Figure ^, one should note that the wind vari- 
ation does not affect the decay law of rik after the jump. 
Actually, the density of kinks follows Eq.(||) before and 
after the wind change. 

Typical fits using Eq. H are drawn in Figure I and 
parameters are reported in Table 1. The upper line of 
Figure ^ correspond to the case of a variable wind, and 
is caracterized by d w 1, while the lower curve is for a 
constant wind and follows d w 2. This difference implies 
a greater stability of the labyrinthic structure, and the 
existence of two modes. 



Constant wind 2.45 x lO"'^ ± 7.6 x 10"'' 1.98 ± 0.03 
Variable wind 9.68 x lO"'^ ± 1 x 10"'' 1.16 ± 0.09 



TABLE I. Parameters a and d of Eq.(y) fitted for both 
constant and variable winds. 

An interesting observation is that if a second wind 
change occurs on the labyrinthic structure, is not af- 
fected. Once the diagonal structure is created, any come 
back to the transveral one is prohibited. The process is ir- 
reversible! However, if the initial topography is composed 
by ripples with a small amount of defects, the landscape 
can evolve to a nearly transversal structure. This behav- 
ior comes from the lack of kinks. Indeed, if is initialy 
small, a few number of labyrinths will be formed. As a 
consequence, the lanscape is less stable. 

Moreover, the formation of labyrinthic structures in- 
duces a kind of "memory effect" . Indee, asymptotic val- 
ues a listed in Table 1 are significantly different if one 
compares constant and variable cases! A difference in 
asymptotic values is a strong result supporting the idea 
that there is a memory of the wind direction on the land- 
scape evolution. After a wind change, the evolution of 
the landscape depends essentialy on the former topogra- 
phy. Even after a long time, the surface always evolves 
in a way depending on its history, i.e on the number on 
wind orientation changes. The question is to know if real 
granular landscapes show this memory effect. This is left 
for future experimental work. 



Since the kink density decreases in a faster way than a 
logarithmic-like law, and in a slower way than an expo- 
nential one, we have assumed a power law decay 

"fcW = « + ^:^ (3) 

where a, &, c and d are fitting parameters. The param- 
eter a represents the asymptotic density of kinks, while 



IV. SUMMARY 

In summary, we have simulated unusual labyrinthic 
landscapes observed in earlier experiments. We have 
investigated the formation and evolution of these land- 
scapes. We have demonstrated that the density of defects 
in a ripple structure is a relevant parameter to character- 
ize the temporal evolution of such structures. Indeed, the 
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number of "defects" present in the landscape decreases 
according to a negative power law of time. If wind orien- 
tation is changed, the power exponent shifts from a value 
2 to the value 1. These exponents do not dependent on 
the occurence of wind change. We have also shown the 
emergence of a memory effect in the asymptotic value of 
the kink density. 
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